Final Presentation

Loading Libraries

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0      ✔ purrr   1.0.0 
## ✔ tibble  3.1.8      ✔ dplyr   1.0.10
## ✔ tidyr   1.2.1      ✔ stringr 1.5.0 
## ✔ readr   2.1.3      ✔ forcats 0.5.2 
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(janitor)

## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(shiny)
library(shinydashboard)

## 
## Attaching package: 'shinydashboard'
## 
## The following object is masked from 'package:graphics':
## 
##     box

library(ggplot2)
library(naniar)
library(lubridate)

## Loading required package: timechange
## 
## Attaching package: 'lubridate'
## 
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(dplyr)
library(RColorBrewer)
library(paletteer)
library(readr)
library(ggmap)

## ℹ Google's Terms of Service: <]8;;https://mapsplatform.google.comhttps://mapsplatform.google.com]8;;>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.

library(devtools)

## Loading required package: usethis

library(here)

## here() starts at C:/GitHub Repositories/GOBY - Edits

library(ggrepel)

Loading & Cleaning Data

visibility <- read_csv("data/CA_visibility_data.csv")

## Rows: 26633 Columns: 22
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (11): Dataset, SiteCode, Date, SiteName, State, ammNO3f_Unit, ammSO4f_Un...
## dbl (11): POC, Percentile, Latitude, Longitude, Elevation, ammNO3f_Val, ammS...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

visibility_tidy <- visibility%>%
  na_if(-999)

Generating Maps (Oonagh)

Step 1: Create a Base Map

1. Find the range of latitude and longitude for our data:

visibility_tidy %>% 
  select(Latitude, Longitude) %>% #note that the lat and long are in degrees decimal
  summary()

##     Latitude       Longitude     
##  Min.   :33.46   Min.   :-124.1  
##  1st Qu.:34.73   1st Qu.:-121.2  
##  Median :37.22   Median :-119.7  
##  Mean   :37.32   Mean   :-119.7  
##  3rd Qu.:38.98   3rd Qu.:-118.1  
##  Max.   :41.71   Max.   :-116.4

2. Set the bounding box range:

lat <- c(33.46, 41.71)
long <- c(-124.1, -116.4)
bbox <- make_bbox(long, lat, f = 0.05)

3. Generate the bounding box

map1 <- get_map(bbox, maptype = "terrain", source = "stamen")

## ℹ Map tiles by Stamen Design, under CC BY 3.0. Data by OpenStreetMap, under ODbL.

ggmap(map1)

Step 2: Adding Points to Base Map

I. Basic map with site locations

ggmap(map1) + 
  geom_point(data = visibility_tidy, aes(Longitude, Latitude), size=1.5) +
  labs(x= "Longitude", y= "Latitude", title="Air Quality Research Site Locations in CA")

  #ggsave("plot.png", width = 5, height = 5)

II.Basic map with site location and by elevation

ggmap(map1) + 
  geom_point(data = visibility_tidy, aes(Longitude, Latitude, colour=Elevation), size=2.5) +
  labs(x= "Longitude", y= "Latitude", title="Air Quality Research Site Locations in CA")

III.Basic map with site location and colour name key - not informative

ggmap(map1) + 
  geom_point(data = visibility_tidy, aes(Longitude, Latitude, colour=SiteName), size=2.5) +
  labs(x= "Longitude", y= "Latitude", title="Air Quality Research Site Locations in CA")

IV.Basic map with only ONE site location - not informative

visibility_yos <- visibility_tidy %>% 
  filter(SiteName=="Yosemite NP")

  ggmap(map1) + 
  geom_point(data = visibility_yos, aes(Longitude, Latitude), size=2.5) +
  labs(x= "Longitude", y= "Latitude", title="Yosemite Air Quality Research Site")

Saving maps

ggsave("plot.png", width = 5, height = 5) #Saves last plot as 5’ x 5’ file named "plot.png" in working directory.
#Matches file type to file extension.

Split date by year:

visibility2 <- visibility_tidy %>% 
  separate(Date, into= c("month", "day", "year"), sep = "/")

visibility3 <- visibility2 %>% 
  group_by(SiteName, year) %>% 
  summarise(across(c(SVR_Val, ammNO3f_Val, ammSO4f_Val, ECf_Val, OMCf_Val, SOILf_Val, Longitude, Latitude), mean, na.rm=T), total=n())

## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.

A) Graph Mean Standard Visual Range over time, faceted by site

visibility3 %>% 
  ggplot(aes(x=year, y=SVR_Val))+
  geom_point(na.rm=TRUE)+
  facet_wrap(~SiteName, ncol=4)+
  theme(axis.text.x = element_text(angle = 60, hjust = 1))+
  labs(title = "Standard Visual Range over Time",
       x = "Year",
       y = "Average Standard Visual Range (km)",) #fill by season?

B) Graph Mean Ammonium Nitrate Concentration Over Time faceted by site

visibility3 %>% 
  ggplot(aes(x=year, y=ammNO3f_Val))+
  geom_point(na.rm=TRUE)+
  facet_wrap(~SiteName, ncol=4)+
  theme(axis.text.x = element_text(angle = 60, hjust = 1))+
  labs(title = "Ammonium Nitrate Concentration over Time",
       x = "Year",
       y = "Ammonium Nitrate Concentration (ug/m^3)",) #fill by season?

C) Graph Mean Ammonium Sulfate Concentration Over Time faceted by site

visibility3 %>% 
  ggplot(aes(x=year, y=ammSO4f_Val))+
  geom_point(na.rm=TRUE)+
  facet_wrap(~SiteName, ncol=4)+
  theme(axis.text.x = element_text(angle = 60, hjust = 1))+
  labs(title = "Ammonium Sulfate Concentration over Time",
       x = "Year",
       y = "Ammonium Sulfate Concentration (ug/m^3)",) #fill by season?

D) Graph Mean Elemental Carbon Concentration Over Time faceted by site

visibility3 %>% 
  ggplot(aes(x=year, y=ECf_Val))+
  geom_point(na.rm=TRUE)+
  facet_wrap(~SiteName, ncol=4)+
  theme(axis.text.x = element_text(angle = 60, hjust = 1))+
  labs(title = "Elemental Carbon Concentration over Time",
       x = "Year",
       y = "Elemental Carbon Concentration (ug/m^3)",) #fill by season?

Notes: Interesting fluctuations for Lassen, Lava Beds, Point Reyes, Yosemite, Sequoia, Trinity. Link to wildfire data?

Generating Wildfire Graphs (Yukari)

Separate Dates for plotting

visibility_mdy <- visibility_tidy%>%
  separate(Date, into = c("month", "day", "year"))
visibility_mdy

## # A tibble: 26,633 × 24
##    Dataset SiteCode   POC month day   year  Percentile SiteName  Latit…¹ Longi…²
##    <chr>   <chr>    <dbl> <chr> <chr> <chr>      <dbl> <chr>       <dbl>   <dbl>
##  1 IMPRHR2 AGTI1        1 1     3     2011          10 Agua Tib…    33.5   -117.
##  2 IMPRHR2 AGTI1        1 1     6     2011          10 Agua Tib…    33.5   -117.
##  3 IMPRHR2 AGTI1        1 1     9     2011          70 Agua Tib…    33.5   -117.
##  4 IMPRHR2 AGTI1        1 1     12    2011          10 Agua Tib…    33.5   -117.
##  5 IMPRHR2 AGTI1        1 1     15    2011          10 Agua Tib…    33.5   -117.
##  6 IMPRHR2 AGTI1        1 1     18    2011          10 Agua Tib…    33.5   -117.
##  7 IMPRHR2 AGTI1        1 1     21    2011          10 Agua Tib…    33.5   -117.
##  8 IMPRHR2 AGTI1        1 1     24    2011          10 Agua Tib…    33.5   -117.
##  9 IMPRHR2 AGTI1        1 1     27    2011          10 Agua Tib…    33.5   -117.
## 10 IMPRHR2 AGTI1        1 1     30    2011          50 Agua Tib…    33.5   -117.
## # … with 26,623 more rows, 14 more variables: Elevation <dbl>, State <chr>,
## #   ammNO3f_Val <dbl>, ammNO3f_Unit <chr>, ammSO4f_Val <dbl>,
## #   ammSO4f_Unit <chr>, ECf_Val <dbl>, ECf_Unit <chr>, OMCf_Val <dbl>,
## #   OMCf_Unit <chr>, SOILf_Val <dbl>, SOILf_Unit <chr>, SVR_Val <dbl>,
## #   SVR_Unit <chr>, and abbreviated variable names ¹​Latitude, ²​Longitude

Make Plots Yosemite National Park - 2016 -> 2017 -> 2018 Empire Fire, Ferguson Fire

visibility_mdy %>% 
  group_by(SiteName, year) %>% 
  filter(SiteName=="Yosemite NP") %>% 
  summarize(mean_ano3_year= mean(ammNO3f_Val, na.rm = T)) %>% 
  ggplot(aes(x=year, y=mean_ano3_year)) +
  geom_line(color="red", 
            alpha = 0.8, 
            group = 1) +
  geom_point() +
  geom_label(label = "Ferguson Fire",
            x = 8,
            y = .59,
            color = "black",
            fill = NA) +
  labs(title = "Mean Ammonium Nitrate in Yosemite National Park",
       x = "Year",
       y = "Ammonium Nitrate (ug/m^3)")

## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.

ggsave("yosemitenp_ammno3.png",
       plot = last_plot())

## Saving 7 x 5 in image

visibility_mdy %>% 
  group_by(SiteName, year) %>% 
  filter(SiteName=="Yosemite NP") %>% 
  summarize(mean_aso4_year= mean(ammSO4f_Val, na.rm = T)) %>% 
  ggplot(aes(x=year, y=mean_aso4_year))+
  geom_line(color="red", 
            alpha = 0.8, 
            group = 1) +
  geom_point() + 
  geom_label(label = "Ferguson Fire",
            x = 8,
            y = .79,
            color = "black",
            fill = NA) +
  labs(title = "Mean Ammonium Sulfate in Yosemite National Park",
       x = "Year",
       y = "Ammonium Sulfate (ug/m^3)")

## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.

ggsave("yosemitenp_ammso4.png",
       plot = last_plot())

## Saving 7 x 5 in image

visibility_mdy %>% 
  group_by(SiteName, year) %>% 
  filter(SiteName=="Yosemite NP") %>% 
  summarize(mean_omc_year= mean(OMCf_Val, na.rm = T)) %>% 
  ggplot(aes(x=year, y=mean_omc_year))+
  geom_line(color="red", 
            alpha = 0.8, 
            group = 1) + 
  geom_point() +
  geom_label(label = "Ferguson Fire",
            x = 8,
            y = 10,
            color = "black",
            fill = NA) +
  labs(title = "Mean OMC in Yosemite National Park",
       x = "Year",
       y = "Organic Mass by Carbon (ug/m^3)")

## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.

ggsave("yosemitenp_omc.png",
       plot = last_plot())

## Saving 7 x 5 in image

Lassen Volcanic National Park - 2020 -> 2021 Dixie Fire

visibility_mdy %>% 
  group_by(SiteName, year) %>% 
  filter(SiteName=="Lassen Volcanic NP") %>% 
  summarize(mean_ano3_year= mean(ammNO3f_Val, na.rm = T)) %>% 
  ggplot(aes(x=year, y=mean_ano3_year)) +
  geom_line(color="red", 
            alpha = 0.8, 
            group = 1) +
  geom_point() +
  geom_label(label = "Dixie Fire",
            x = 11,
            y = .22,
            color = "black",
            fill = NA) +
  labs(title = "Mean Ammonium Nitrate in Lassen Volcanic National Park",
       x = "Year",
       y = "Ammonium Nitrate (ug/m^3)")

## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.

ggsave("lassenvolcanicnp_ammno3.png",
       plot = last_plot(),
       width = 7,
       height = 5)

visibility_mdy %>% 
  group_by(SiteName, year) %>% 
  filter(SiteName=="Lassen Volcanic NP") %>% 
  summarize(mean_aso4_year= mean(ammSO4f_Val, na.rm = T)) %>% 
  ggplot(aes(x=year, y=mean_aso4_year))+
  geom_line(color="red", 
            alpha = 0.8, 
            group = 1) +
  geom_point() + 
  geom_label(label = "Dixie Fire",
            x = 11,
            y = .465,
            color = "black",
            fill = NA) +
  labs(title = "Mean Ammonium Sulfate in Lassen Volcanic National Park",
       x = "Year",
       y = "Ammonium Sulfate (ug/m^3)")

## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.

ggsave("lassenvolcanicnp_ammso4.png",
       plot = last_plot(),
       width = 7,
       height = 5)

visibility_mdy %>% 
  group_by(SiteName, year) %>% 
  filter(SiteName=="Lassen Volcanic NP") %>% 
  summarize(mean_omc_year= mean(OMCf_Val, na.rm = T)) %>% 
  ggplot(aes(x=year, y=mean_omc_year))+
  geom_line(color="red", 
            alpha = 0.8, 
            group = 1) + 
  geom_point() +
  geom_label(label = "Dixie Fire",
            x = 11,
            y = 5.6,
            color = "black",
            fill = NA) +
  labs(title = "Mean OMC in Lassen Volcanic National Park",
       x = "Year",
       y = "Organic Mass by Carbon (ug/m^3)")

## `summarise()` has grouped output by 'SiteName'. You can override using the
## `.groups` argument.

ggsave("lassenvolcanicnp_omc.png",
       plot = last_plot(),
       width = 7,
       height = 5)

Shiny App: Mean Standard Visual Range per Year by Site (Bao)

improve_visibility <- visibility %>% 
  clean_names() %>% 
  na_if("-999") %>% 
  select(date, site_name, state,svr_val) %>% 
  separate(date, into = c("month", "day", "year"))

ui <- dashboardPage(
  dashboardHeader(title = tags$div(style = "font-size: 30px;",
    "Visual Range CA")),
  dashboardSidebar(disable = TRUE),
  
  dashboardBody(
      tags$head(
      tags$style(HTML(".box{border-radius: 20px; box-shadow: 5px 5px 10px #888888;}"))),
      
  fluidRow( class="text center",
  box(title= "Plot Options",width = 2,
      radioButtons("x", "Choose The Site", choices = c("Agua Tibia", "Bliss SP (TRPA)", "Bliss SP (TRPA) (RHTS)", "Death Valley NP", "Dome Lands Wilderness", "Fresno", "Hoover", "Joshua Tree NP", "Kaiser", "Lake Tahoe Community College", "Lassen Volcanic NP", "Lava Beds NM", "Owens Valley", "Pinnacles NM", "Point Reyes National Seashore", "Redwood NP", "San Gabriel", "San Gorgonio Wilderness", "San Rafael", "Sequoia NP", "Trinity", "Wrightwood", "Yosemite NP")
)
),  #Close box     

box(title = h2("Standard Visual Range Over the Year", align="center"), width = 6, align="center",solidHeader = T,
plotOutput("plot", width = "550px", height = "525px")),

box(title="", width = 4, height = 650,
imageOutput("image"), style ="text-align: center;"
),
) #close fluidrow
) #close dashboardbody
) #close dashboard page

server <- function(input, output, session) {
  session$onSessionEnded(stopApp)
  output$plot <- renderPlot({
  improve_visibility %>% 
  group_by(site_name, year) %>% 
  summarize("Mean_Standard_Visual_Range_km"= mean(svr_val, na.rm = T)) %>% 
  filter(site_name==input$x) %>% 
  ggplot(aes(x=year, y=Mean_Standard_Visual_Range_km,group=2, color=Mean_Standard_Visual_Range_km))+
  geom_line(size=1.5,alpha=0.8)+
  geom_point(size=5)+
      theme(axis.title.x = element_text(size = 23, margin = margin(t = 15)))+
      theme(axis.title.y = element_text(size = 23,  margin = margin(r = 20)))+
      theme(axis.text = element_text(size = 14))+
      labs(y = "Average Standard Visual Range (km)")+
      labs(x = "Year")+
      scale_color_gradient(low="blue", high="red")+
      guides(color = FALSE)
  })
  output$image <- renderImage({
    list(src = "map.png", height = "570px", width="450px")
  }, deleteFile = FALSE)
    
  }

shinyApp(ui, server, options = list(launch.browser = TRUE))

Shiny App: Seasonal Metrics per Year by Site (Gabe)

Designating Seasons Warning: This code block is not time efficient (O(n)), expect a long run time.

seasons <- c()

for (date in visibility_tidy$Date){
  date <- mdy(date)
  # Winter Dec 20 - Mar 19 
  if ((month(date) == 12 & day(date) >= 20) | month(date) %in% c(1, 2) | (month(date) == 3 & day(date) <= 19)){
    seasons <- c(seasons, "winter")
  }
  # Spring Mar 20 - June 20
  else if ((month(date) == 3 & day(date) >= 20) | month(date) %in% c(4, 5) | (month(date) == 6 & day(date) <= 20)){
    seasons <- c(seasons, "spring")
  }
  # Summer June 21 - Sept 22
  else if ((month(date) == 6 & day(date) >= 21) | month(date) %in% c(7, 8) | (month(date) == 9 & day(date) <= 22)){
    seasons <- c(seasons, "summer")
  }
  # Autumn Sept 23 - Dec 19
  else{
    seasons <- c(seasons, "autumn")
  }
}

visibility_tidy$season = seasons

Separate Dates for plotting

visibility_tidy_mdy <- visibility_tidy%>%
  separate(Date, into = c("month", "day", "year"))
head(visibility_tidy_mdy)

## # A tibble: 6 × 25
##   Dataset SiteCode   POC month day   year  Percentile SiteName   Latit…¹ Longi…²
##   <chr>   <chr>    <dbl> <chr> <chr> <chr>      <dbl> <chr>        <dbl>   <dbl>
## 1 IMPRHR2 AGTI1        1 1     3     2011          10 Agua Tibia    33.5   -117.
## 2 IMPRHR2 AGTI1        1 1     6     2011          10 Agua Tibia    33.5   -117.
## 3 IMPRHR2 AGTI1        1 1     9     2011          70 Agua Tibia    33.5   -117.
## 4 IMPRHR2 AGTI1        1 1     12    2011          10 Agua Tibia    33.5   -117.
## 5 IMPRHR2 AGTI1        1 1     15    2011          10 Agua Tibia    33.5   -117.
## 6 IMPRHR2 AGTI1        1 1     18    2011          10 Agua Tibia    33.5   -117.
## # … with 15 more variables: Elevation <dbl>, State <chr>, ammNO3f_Val <dbl>,
## #   ammNO3f_Unit <chr>, ammSO4f_Val <dbl>, ammSO4f_Unit <chr>, ECf_Val <dbl>,
## #   ECf_Unit <chr>, OMCf_Val <dbl>, OMCf_Unit <chr>, SOILf_Val <dbl>,
## #   SOILf_Unit <chr>, SVR_Val <dbl>, SVR_Unit <chr>, season <chr>, and
## #   abbreviated variable names ¹​Latitude, ²​Longitude

Explore the functionality of separating by season Finding the mean data per year per season.

mean_by_season <- visibility_tidy_mdy%>%
  mutate(SiteName=paste(SiteName, year, sep="_"))%>%
  mutate(SiteName=paste(SiteName, season, sep=" "))%>%
  group_by(SiteName)%>%
  summarise("mean_ammNO3_ppb"=mean(ammNO3f_Val, na.rm=T),
            "mean_ammSO4_ppb"=mean(ammSO4f_Val, na.rm=T),
            "mean_EC_ppb"=mean(ECf_Val, na.rm=T),
            "mean_OMC_ppb"=mean(OMCf_Val, na.rm=T),
            "mean_SOIL_ppb"=mean(SOILf_Val, na.rm=T),
            "mean_SVR_km"=mean(SVR_Val, na.rm=T))%>%
  as.data.frame()%>%
  separate(SiteName, into = c("SiteName", "year"), sep="_")%>%
  separate(year, into = c("year", "season"), sep=" ")
head(mean_by_season)

##     SiteName year season mean_ammNO3_ppb mean_ammSO4_ppb mean_EC_ppb
## 1 Agua Tibia 2011 autumn       0.7042493        1.236351   0.2399840
## 2 Agua Tibia 2011 spring       1.3284576        1.929723   0.1828920
## 3 Agua Tibia 2011 summer       0.9324993        2.694624   0.2530929
## 4 Agua Tibia 2011 winter       0.6691160        0.510433   0.1339632
## 5 Agua Tibia 2012 autumn       0.6549057        1.049741   0.1992536
## 6 Agua Tibia 2012 spring       2.0285125        2.128372   0.2105500
##   mean_OMC_ppb mean_SOIL_ppb mean_SVR_km
## 1    1.5392304     0.4572863   137.75418
## 2    1.7740872     0.8754060    93.61895
## 3    2.2604336     0.5468550    84.92630
## 4    0.8716737     0.2251410   191.02952
## 5    1.2527293     0.4292511   141.15836
## 6    1.7745840     0.7293125    82.39903

Creating Plots Here is our reference graph:

mean_by_season%>%
  filter(SiteName=="Fresno")%>%
  ggplot(aes(year, mean_SVR_km, fill=season))+
  geom_col(position = "dodge")

Here is the shiny app:

ui <- dashboardPage(
  dashboardHeader(title="Seasonal Metrics"),
  dashboardSidebar(disable = T),
  dashboardBody(
    
    tags$head(
      tags$style(HTML(".box{border-radius: 20px; box-shadow: 5px 5px 10px #888888;}"))),
    
    fluidRow( 
      class="text center",
      box(
        selectInput("site", "Select Site Name", choices = unique(visibility_tidy_mdy$SiteName))
        ),
      box(
        selectInput("y", "Visibility Metric", 
              choices = c("Standard Visual Range"="mean_SVR_km", "Ammonium Sulfate" = "mean_ammSO4_ppb", 
                          "Ammonium Nitrate" = "mean_ammNO3_ppb", "Elemental Cabon" = "mean_EC_ppb", 
                          "Organic Mass Cabon" = "mean_OMC_ppb", "Dust" = "mean_SOIL_ppb"),
              selected = "mean_SVR_km")
      ),
      box(
         plotOutput("plot", width = "1400px", height = "400px"), width = 12
      )
    )
  )
)

server <- function(input, output, session) {
  session$onSessionEnded(stopApp)
  output$plot <- renderPlot({
    mean_by_season%>%
      filter(SiteName == input$site)%>%
      ggplot(aes_string("year", input$y, fill="season"))+
      scale_fill_manual(values=c("#00CCE6", "#9BEB01", "#E6E500", "#EB7500"))+
      geom_col(position = "dodge", color="black")+
      theme(axis.title = element_text(size = 16),
            legend.text = element_text(size = 14),
            legend.title = element_text(size = 16),
            legend.background = element_rect(fill = "#F1F1F1"))+
      labs(x="Year", 
           y=paste(str_split_1(input$y, "_")[2], str_split_1(input$y, "_")[3], sep=" "),
           fill="Season")
  })
}

shinyApp(ui, server, options = list(launch.browser = TRUE))